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Abstract. We present the Parareal-CG algorithm for time-dependent differential equations in this work. 
The algorithm is a parallel in time iteration algorithm utilizes Chebyshev-Gauss spectral collocation method 
for fine propagator F and backward Euler method for coarse propagator G. As far as we know, this is the 
first time that the spectral method used as the F propagator of the parareal algorithm. By constructing 
the stable function of the Chebyshev-Gauss spectral collocation method for the symmetric positive definite 
(SPD) problem, we find out that the Parareal-CG algorithm and the Parareal-TR algorithm, whose F prop- 
agator is chosen to be a trapezoidal ruler, converge similarly, i.e., the Parareal-CG algorithm converge as 
fast as Parareal-Euler algorithm with sufficient Chebyhsev-Gauss points in every coarse grid. Numerical ex- 
amples including ordinary differential equations and time-dependent partial differential equations are given 


to illustrate the high efficiency and accuracy of the proposed algorithm. 
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1 Introduction 


We are considering utilizing the parareal algorithm for the initial-value problems in the following 


du 
T = f(t,u), te (0, T], (1.1) 
u(0) = uo 


where f : (0,7) xR™ — R”, and up € R™. Parareal is a well-studied parallel in time algorithm developed by 
Lions, Maday, and Turinici in 2001 [27]. The algorithm obtains the solution in a limited number of predictor- 
corrector iterations utilizing random initial values at each temporal subinterval, stopping when a tolerance 


is reached. The global error produced by this iterative method is equivalent to that obtained by the serial 
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of the fine propagator. Due to the advantages of parareal algorithm including but not limited to efficiency 
and convergency, many relevant methods have emerged in recent years [14, 27, 31, 43]. Together with these 
methods, the parareal algorithm has been identified in various fields of research, including optimal control 
problems [13, 30, 36, 38], wave equations [11, 20], stochastic differential equations [5, 22, 46], Hamiltonian 
systems [9, 15], incompressible flows [8, 10], heat equations [23, 25, 40], algebraic equations [16, 24], molecular 
dynamics [2, 26], and other partial differential equations [4, 28, 29]. 

The parareal algorithm combines two numerical methods which are the coarse propagator G and fine 
propagator F, associated with the large time step AT and small time step At, respectively. The ratio 
J = AT/At is assumed to be greater than 1. The G propagator is often chosen to be backward Euler method 
which is inexpensive and strongly stable so that it is available for the large time step AT computations. 
Additionally, numerical methods based on Taylor’s expansion and quadrature formula which is much cheaper 
are frequently chosen as F propagators, and the convergence of various F propagators has been thoroughly 


investigated for the symmetric positive definite (SPD) problem. 
u' (t) + Au(t) = g(t), A€ R™*™ symmetric positive definite (SPD). (1.2) 


Gander and Vandewalle [17] proposed the convergence theorem of parareal algorithm and illustrated numer- 
ically that Parareal-Euler algorithm using for F the backward-Euler method converges rapidly. Based on 
which, Mathew, Sarkis and Schaerer [30] proved theoretically that the parareal-Euler algorithm converged 
robustly and the convergence factor is 0.298 for J > 2. Wu [35] demonstrated that the robust convergence 
also holds for F propagator chosen to be second-order diagonal implicit Runge-Kutta (DIRK2) method 
and TR/BDF2 method (ode23tb solver for ODEs in Matlab). For J > 2, the convergence factor of the 
Parareal-2sDIRK and Parareal-TR/BDF2 are 0.316 and 0.333 respectively. Wu and Zhou [42] also showed 
the analysis F propagator for the third-order diagonal implicit Runge-Kutta method with a convergence 
factor 0.333 for J > 4; they also proved that for trapezoidal formula and fourth-order Gauss-Runge-Kutta 
integrator chosen to be F propagator, there exists a Jž in depends on both spectral radius of A and the step 
size At which make the convergence factor be 0.333 for J > J* 


žin: Recently, Yang, Yuan and Zhou [45] gave 
a more general result stating that if the F propagator is strongly stable single step integrators, there must 


exist a positive J*,,, (independent of step sizes AT, At, terminal time T, problem data uo and f, as well as 


min 
the spectral radius of A) such that the parareal algorithm converges linearly with convergence factor close 
to 0.3 for all J > Jýin- 

In this work, we choose Chebysheyv-Gauss spectral collocation method [44] for fine propagator and present 
Parareal-CG algorithm. The Chebyshev-Gauss spectral collocation method is an overall iteration method 
with M Chebyshev-Gauss points in the computation interval, so it combines the advantages of spectral 
accuracy and computational efficiency. Also, it is unnecessary to solve implicit equations for every fine 
time step At as classical methods via Taylor’s expansion mentioned in [17, 35, 42, 45]. We would briefly 
introduce the spectral collocation methods for ODEs. Clenshaw and Norton first presented the Chebyshev- 
Picard method for solving nonlinear ordinary differential equations in 1963, in which the nonlinear term was 
approximated by Chebyshev series, so that the method was collocated at Chebyshev-Gauss-Lobatoo points 


and implemented by Picard iteration. Later, a matrix-vector form of the method was introduced by Feagin 


and Nacozy [12], greatly increases the computation efficiency. Yang and Wang [44] proposed Chebyshev- 
Gauss spectral collocation method via Chebyshev-Gauss points for ODEs in a single interval and analyzed 
the convergence by the hp version. The method demonstrated significant advantages in astrodynamics 
simulations [1, 3, 37, 41] because of its spectral accuracy and computational efficiency. Additionally, some 
relevant spectral collocation methods for solving ODEs are also proposed in recent years [18, 19, 39]. 

We present the matrix-vector form of Chebyshev-Gauss spectral collocation method by the approach in 
[12]. Based on which we construct its stability function and illustrate the convergence of the Parareal-CG 
method numerically. We find out there exists a Mžin depends on spectral radius of A and the step size AT 
which makes the convergence factor be 0.333 for M > Mă 


min’ 


The proposed method has the same convergence 
as Parareal-TR and Parareal-Gauss4 methods. 

The rest of the paper is organized as follows, in section 2 we recall the Chebyshev-Gauss spectral collo- 
cation method and its convergence theorem. Parareal-CG algorithm is proposed in 3. We provide stability 
function of Chebyshev-Gauss spectral collocation method and observe the convergence of Parareal-CG al- 
gorithm in 4. Several numerical experiments are carried out in Section 5 to demonstrate the high accuracy 


and convergency of the proposed method. We finally give some conclusions in Section 6. 


2 Revisit of Chebyshev-Gauss Spectral Collocation Method 


In this section, we revisit the Chebyshev-Gauss spectral collocation method proposed by Yang and Wang 
[44] in a general interval [a,b], (b > a > 0) with the initial condition u(a) = ua. 

Let Tı(T) = cos(larccos(T)) be the standard Chebyshev polynomials of degree l, (l = 0,1,---) with 
7 € [-1,1], then by using the affine transformation we can define the shifted Chebyshev polynomials 


fw =n (E21), te a,b], §=0,1,--- (2.1) 


According to the definition, one can obtain the following shifted Chebyshev derivative relationship directly 


~ 2 
T(t) = 0, Ti (t) = T, 
i ~~“ (2.2) 
<y Ti (t D, 123 
yi me) l 1 1-1(t) b— 1( ); 
Let Tm denote the standard Chebyshev-Gauss (CG) points in (—1, 1), 
(2m + 1)r 
eee = ia 2. 
Tin cos ioga? ™ 0,1,- ,M, (2.3) 
the corresponding shifted Chebyshev-Gauss (CG) points tm has the form 
b— b 
tm = rn — y m=0,1,:-,M, (2.4) 


which are the zeros of Tija A): 


Denote Pm+ı(a,b) be a set of polynomials of degree at most M + 1 in (a,b), the Chebyshev-Gauss 
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spectral collocation method is to seek um(t) € Pm+1ı(a, b) defined by 


M-+1 a 
m=0 
such that 4 
war m =7 ms m)); = ,1,---,M, 
gí ) uf (t um (t )) m=0 (2.6) 


um(a) = Ua, 


where Zy f (t, u(t)) : C(a, b) + Par(a, b) is the Chebyshev interpolation of f (t, u(t)) defined by 
M 
Im fit, um (t)) = 5 ÍmTm(t), (2.7) 
m=0 


the coefficients { Ta M_, are determined by the forward discrete Chebyshev transform 


M 
fa = Sag Mma) Flt) (2.8) 


where co = 2, Cm = 1 for m > 1. Then by the shifted Chebyshev derivative relationship (2.2), we can derive 
the coefficients {ûm} 1a in (2.5) by 


(b-a) . (b-a) 


UM+1 = IM UM = 4M fm-1, 

. b—a A z 

Um = j emai din—1 = fmi); 1 < m < M- 1; (2.9) 
M+1 

to = Ua — 5 (—1)*ti,. 
k=1 


Yang and Wang [44] presented the error estimate for the hp-version of the single interval Chebyshev-Gauss 
spectral collocation method, before introduce the theorem, we first present some notations used throughout 


the error estimate 


e H? (a,b), (r > 0) denotes the weighted Sobolev space with certain weight function w = (t — a)(b-— t) 
in (a,b), especially, H°(a,b) = L? (a, b). 


e ||vl|. denotes the norm of the space L?_,,, (a,b). 
The error estimate theorem is stated as follows. 


Theorem 2.1 Assume that f(t, z) fulfills the Lipschitz condition, 
Frnt) — Fzt) < Lla = zl, L>0, (2.10) 


and 0 < L(b—a) < 8 < } (B is a certain constant) holds. Then for any u € H” (a,b) with integers 
Ww 2 


2<r<N-+1, we have 


b-a a d” 5 
lu — umllt2(a,b) < 5 lu —uarl, < Coo- ayer f w" 2 (t) (Fu) dt, 


a 


b . 2 
ju) — uar(0)] < Ca(0— aya f uri (Tate) a 


—u 
ie dt” 
where Cg is a positive constant depending only on p. 


In actual computation, we use Picard iteration procedure to compute the coefficients {ti }44_,. The 
approach is simple to implement, especially for the complex nonlinear problems. The p-th (p = 1,2,---) 


Picard iteration form of (2.6) is 


d 
TUM (tm) = Tuf (tm ths (tm)) (2.12) 


if the equation satisfies the convergence condition in Theorem 2.1, the iteration solution u4,(t) will converge 
to the numerical solution um (t) with large enough p, and the convergence is of order one, that is to say there 


exists a constant 0 < Cp < 1 such that |lep+1|]oo < Cpllep||.o with the definition ep := up! (t) —u™ (t). 


Algorithm 1 Chebyshev-Gauss Spectral Collocation Algorithm 


Input: Provide the initial guess of {u (tm) }_9, the tolerance e. 
For p = 0,1,- 
Step 1. Evaluate the values of {f (tm, up! (tm))}u=o: 


Step 2. Compute the coefficients {fin tro by (2.8). 
Step 3. Compute the coe Cne {0 m=0 by (2.9). 
Step 4. Update the data of {uf} (tm }m=0 by (2-5). 
Step 5. If the iteration error satisfies the stopping criterion, 


1 
luht (t) — uir Ollo < €, 


terminate the iteration; otherwise go back to Step 1. 
Step 6. Compute ult (b) _ yi M41 pp 


m=0 “mn: 


We designate Fca as the numerical propagator defined by the Chebyshev-Gauss spectral collocation 
method at the conclusion of the section. The numerical output of Algorithm 1 with M points in the interval 


[a, b] is represented as Foq(a, ua, M,b — a). That is, 


uh (b) = Foa(a, ta, M,b— a). (2.13) 


3 Parareal-CG algorithm 


The parareal algorithm introduced by Gander and Vandewalle [17] is revisited in this section, using the 
Chebyshev-Gauss spectral collocation method as the fine propagator. First, we divide the whole time interval 
(0, T] uniformly by 0 = To < T; <---< Ty =T and define AT = T/N. Second, we divide each (Tn, Tn+1) 
by M(> 2) shifted Chebyshev-Gauss points (2.4). Then, the low-order and inexpensive numerical method G 


propagator is applied to the coarse time grids, while the Chebyshev-Gauss spectral collocation method Fee 


propagator having spectral accuracy is utilized in the fine grids. The time-sequential and time-parallel parts 
of the algorithm are denoted by the symbols © and 9, respectively. The following is the parareal method 
employing Fag as the fine propagator. 


Algorithm 2 Parareal-CG Algorithm 


© Initialization: Compute sequentially u?,, = G(T,,u?, AT) with u8 = uo, n=0,1,---,N—-1; 
For k = 0,1,- 
® Step 1. On each subinterval [Tn, Tn+41], compute ŭn+ı = Fca (Tn, uk, M, AT). 
© Step 2. Perform sequential corrections 
uppi = G(In Unt, AT) + ny1 — G (Th, Uy, AT), (3.1) 
where ukt! = ug, n=0,1,---,N—1; 
© Step 3. If {u**+1}4_, satisfies the stopping criterion, terminate the iteration and output {u8 t1}N_]; 


n=1 


otherwise go back to Step 1. 


Given that the implicit Euler method is L-stable and feasible for high coarse time step sizes AT forced 
by the parareal algorithm, using it as the coarse propagator makes sense. Other reliable implicit numerical 
methods, such implicit Runge-Kutta methods, are also available for use as the G propagator, although they 
are all significantly more costly than the implicit Euler method. In this paper, we apply the Chebyshev- 
Gauss spectral collocation method to the fine propagator, and the compact form of the Parareal-CG method 
may be expressed as 


uk ti =G(Tn, unt’, AT) + Foc (Tn, uk, M, AT) — G (Th, U 


n? 


E ATY: (3.2) 


n? 

For comparison, we also take into account other high-order and expensive numerical methods like trape- 
zoidal formula. Accordingly, each interval |In, Tn+1] should be divided into J (> 2) small time-intervals 
IT AT 


iT i41], j =0,1,--- ,J—1. We assume the intervals are of uniform size, therefore, At = “+. After 
ntg’ n+ > J 


that, the parareal algorithm can be derived by 


Algorithm 3 Parareal Algorithm 


O Initialization: Compute sequentially u? = G(Tn,u?, AT) with uj = uo, n =0,1,---,N—-1; 

For k = 0,1,- 

® Step 1. On each subinterval (Tr, In41], compute thy itt = F(T pis tng ds Ar) with initial value 
ün = uk, where Ta} 3 =T, + SF and j =0,1,---,J-1; 

© Step 2. Perform sequential corrections 


Un = G(Th, upr, AT) + Un+1 a G(Tn, už, AT), (3.3) 
where ujt? = up, n =0,1,--- , N — 1; 
© Step 3. If {uš+t!}N_] satisfies the stopping criterion, terminate the iteration and output {uk*+1}4_; 


otherwise go back to Step 1. 


The compact form of the parareal algorithm can be derived by 


Unyi = G(Tns tin, AT) + F” (Ta, up, At) — G (Tn, up, AT), oe) 


where F” (Ta, uE, At) stands for the result of running J steps of the fine propagator F with initial value u% 


and the small step-size At. 


4 Convergence Analysis 


Based on the parareal convergence analysis given by Gander and Vandewalle [17], the convergence of 
Parareal-CG method is analysed by constructing the stability function of Chebyshev-Gauss spectral col- 


location method in this section. 


4.1 Stability function of Chevyshev-Gauss Spectral method 


In order to obtain the convergence factor K(z, M) of Parareal-CG algorithm in Algorithm 2, in this subsection 


we present the stability function of Chebyshev-Gauss spectral collocation method Rc«q(z, M). 


Lemma 4.1 For given M and z := AAT, the stability function Rcg(z,M) of Chebyshev-Gauss spectral 
collocation method in [Tn,Tn4i], n = 1,2,--- , N) is 


Rog(z, M) = T (Ie — zCa(Ih + 2T1Coa) 'T,)E, (4.1) 


where Tz and Ca are the coefficient matrices defined in (4.7) and (4.9), respectively; Iı and Iz are two identi- 
ty matrices of the size (M +1) and (M +2), respectively; T = [1,1,--- ,1] and E = [1,0,--- ,0]" 


1x(M+2) (M+4+2)x1 


are two constant vectors. 


Proof. We first introduce the matrix-vector form of the Chebyshev-Gauss spectral collocation method to 
be able to express Reg explicitly. Denote the evaluations of the approximated polynomial for the p-th 
(p = 1,2,---) iteration u4,(t) at the shifted Chebyshev-Gauss points {tm} _o, (n =0,1,--- , N — 1) by 


Up = [uM (to), uM (ti) uM (tu) iuri (4.2) 


Equation (2.9) indicates that the vector up can be expressed as a consequence of the fact Ti(tm) = Ti(t™m), 
(b= 1,2,++*; m =0,1,--- ,M) 


Up = Titi, (4.3) 
where T; is a coefficient matrix defined by 
(M+1) x (M42) 
To(to) Tilto) + Tar+i(t0) 
Ti T) see Ty 
ne olmi) Tifa) M+1(T1) (4.4) 
To(ta) Tilt)  Tm+ı(T™m) 
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Denote f (Gast m) = f®, and suppose the initial value be u(T,) = ur,, then the coefficient vector &? 


can be expressed by 


à? > (ti, a, a ilonsa 
M+1 fp—1 fp—1 
sisa AT ye goed pp—1 AT(fm-2— fm ) AT gp-1 AT 
= —1)"""a?,, af — T M 
UTn + Dut ) Um, 4 ( fo 2 ), ’ 4(M — 1) "AM MW 4(M +1) M 
R:(M+2)x(M+2) S:(M+2)x(M+1) 
UTn 1 2 -4 S2 83 SM 
eo 
0 
0 1 2 0 —1 0 0 
apai 
1 
0 i 0 1 0 -1 0 
p—1 
akela ? 
4 ai 
3 
0 wo 0 0 1 0 -1 
0 + 0 0 0 1 0 
yai 
M 
0 wa} {09 0 0 0 1 
AT * 
= Uo + g ESP., 
where R and S are the coefficient matrices and the first line of S' satisfies 
1 1 
= 1)” , m=2,3,---,M. 
m= (1) (= 1 m- i) í 


Uo and fP can be defined by 


Uo = [ur,,0,0,--- „007 


(M+2)x1? 
fP = [ pp—1 fp-1 ae 
O o41 3 JM lim+1)xi’ 


(4.6) 
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For simplicity, we set the values of the function f(t, u?~!(t)) on the Chebyshev-Gauss points by the notation 
=F (tm u1 (tm)); m =0,1,--- , M, then (2.8) can be expressed as 


apa (fo Tolo) + fP Talr) +- + fhr Tolrm)) 


aa nA n= aah 


uaa T(t) tI Tu (11) +++ + fir Tut) 


V:(M+1)x(M+1) T2:(M+1)x(M+1) (4.7) 
HHI Tolto) Tolna) = Tolrm) | | 
7 ua Tilt) Tifa) = Bila) p= 
ma| (Tuo) Tul) © Tulrm)| |I 
— Vif? . 
where V and T are coefficient matrices and the vector f?—1 is defined by 
f? = [Fat Go) Ft uir Ges , Ftu, way (te) re gaya’ (4.8) 


To this end, given the initial guess ug, the p-th matrix-vector iteration version of the Chebyshev-Gauss 


spectral collocation method for solving (2.6) is as follows. 


1 
Ca = 7RSVT2, 
ûP = Uo — ATCaf?—}, (4.9) 


uP = Ti ûP. 


We develop the matrix-vector form of the method in order to acquire the stability function on the one hand, 
and to considerably improve computing efficiency so that it may be employed in our numerical experiments 
on the other. The matrix-vector form of the method improves the efficiency of computation sufficiently. 
Following, we establish the linear stability function for the Chebyshev-Gauss spectral collocation method, 


for which we initially derive the special form of the method for solving the linear equation, 


uP = T,a? = Ti (Uo — ATCa f?—') 


(4.10) 
= Tı Uo = AATT,CyuP—?. 


By assuming u* = [u*(to),u*(t1),--- ,u*(taz)]' be the convergence solution after a sufficient number of 


iteration, we may rewrite (4.10) for given z := AAT as 


uu = Tı Uo — zT Cau”, 


w (4.11) 
= (Iı +2TıCa) TıUo, 


where Jı is an identity matrix of size (M + 1), then we can derive the corresponding coefficient vector 


až 


a* = [uj,ud,---,u4,4,]" directly by (4.9) 
û* = Uo — ATCa(àu*) = Uo — zCa (I1 + 2T1Ca) T1 U0. (4.12) 


Finally, the stability function of Chebyshev-Gauss spectral collocation method Rce(z, M) could be acquired 


by the compression relations 


*(Tan) _ Ta” 


Rcalz, M) = Ë =T(h — zCa(lı + 211a) Tı )E, (4.13) 
UT, UT, 
where Iz is an (M-+2) dimensional identity matrix, T and E1 are two vectors defined by T = [1,: -+ ,1, 1, cary) 
_ T 
and E=[1,0,--- 01,4... M 


We can obtain the stability functions Reg, particularly when M = 0 and M = 1. 


2—z z? — 8z +16 
R Iis eee =, 
Ip oeU apea 


Rca(z,0) = (4.14) 


For M = 0,1,2,4,20, we displayed |Rce(z, M)| as the functions of zmax := AT Amax in Figure 4.1, which 


shows that for various M we have 


lim [Ree(2)| = 1. (4.15) 


As a result, we can obtain the contraction factor of Parareal-CG method and derive the convergence 


analysis. 


4.2 Convergence analysis of Parareal-CG method 


The following conclusion regarding the convergence of Parareal-CG Algorithm on sufficiently long time 
intervals can be deduced directly from the analysis provided by Gander and Vandewalle [17] for the linear 
system of ODEs (1.2) with A € R™*™ and symmetric positive definite matrix (which therefore can be 


diagonalized and all the eigenvalues are positive real numbers). 


Theorem 4.1 ({17]) Let Fog be the Chebyshev-Gauss spectral collocation propagator with the stability func- 
tion Roa(z) (4.13), and let o(A) = {A1,--: , Àm} be the set of eigenvalues of the matriz A in (1.2). Then, 


the errors {e£} of the parareal-CG algorithm at k-th iteration using the backward-Euler method as the coarse 


10 
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10° 
Zmax(= AT Amax) 


Figure 4.1: The contraction factor [Rca] as a function of zmax, for different M. 


propagator satisfy 
sup ||Veflloo < p* sup ||Veplloo, p= ymax K(ATA, M), (4.16) 
n n €o 


where p is the convergence factor of the parareal algorithm, k > 1 is the iteration indez, and V € R™*™ 
consists of the eigenvectors of A (i.e.,.V~'AV = diag(Ay,:--,Am)). The argument Koa, which is the 
convergence factor corresponding to a single eigenvalue (or in short ’contraction factor” hereafter), is defined 


by 


[Rea(z, m) sa 


Kealz, J (4.17) 


=] 


Alh 
Trz 

The convergence theorem 4.1 states that the behavior of the Parareal-CG algorithm over a long time 
interval is determined by the convergence factor p(M) defined by p(M) := maxzej0,zmax] (2, M) where 
Zmax = AmaxAT and the quantity Amax denotes the maximal eigenvalue (or an upper bound) of the coefficient 
matrix A in (1.2). We can infer from the equation (4.16) that the smaller p(M) is, the faster the algorithm 
converges. To maintain the convergence rate, the convergence factor p(M) prefers to be around z. 

The convergence factor of Parareal-CG algorithm Kce(z, M) as functions of zmax = AT Amax for various 
M are shown in Figure 4.2. We may infer from the behavior of the Kce(z, M) which is similar to the 


convergence factor of Parareal-TR and Parareal-Gauss 4 algorithms discussed in [42] that 
lim Koc (2) = 0, jim Kea(z) = 1. (4.18) 


It implies that the Parareal-CG algorithm cannot keep the convergence rate for arbitrarily large z, however, 
by solving an appropriately designed M that is large enough, it is still possible to get a feasible parareal 
solver. In further detail, as M grows, the convergence factor p(M) := max,¢j9,z,,,,] oa(z,M) tends to be 
around ¿£ for some given z. 

The Convergence analysis can be derived by the convergence theorem 4.1 of Parareal-CG algorithm and 


the stability function (4.1) of Chebyshev-Gauss spectral collocation method. 


11 
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0 i 
10° 10° 10° 104 
Zmax(= ATAmax) 


Figure 4.2: The contraction factor Kcg as a function of Zmax, for different M. 


Theorem 4.2 Given a fixed constant zmax = ATXmax > 0, there exists some positive integer Mžiņ such 
that. 
1 
max Koc(z,M)<-=, if M > Mžin (4.19) 
z€[0,2max] 3 


where Koa is the contraction factors of the Parareal-CG algorithm defined by (4.17) with the stability func- 
tions of Chebyshev-Gauss spectral collocation method Rca given in (4.13). 


The Parareal-CG algorithm’s lower bound, M;,,,,, is provided by 


0 if z < 2, 
Mnin = {1 if 2.< 25 zí, (4.20) 


Méq otherwise. 


where Méq > 1 depends on Zmax and is the minimum positive integer which satisfies the following equation 


3 max 
[Roe (2max,M)| < 3 aie (4.21) 


(Fama) 
Moreover, the quantity zğ is the unique positive root of Kcea(z, 0) = 3 and zj is the maximum positive root 
of Koa(z, 1) = 4, that is to say, 23 = 1 and zł = 8 + 6v2. 

Proof. When M = 0, we can derive the Roa(z,0) and Kee(z, 0) by 
2 


= 2 FA 
=. = — . 4,22 
Rea 242° Kea(z, 0) Fae (z > 0) ( ) 


12 
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It is obvious that Koq(z, 0) is an increasing function with respect to z > 0 and Keaq(z,0) = z has the unique 


root z0* = 1. that is to say, 


1 
3? z € (0, 25] 
Paes Kool, 0) < Zmax i (4.23) 
APEA eas. ZE (25; Zmaxl: 
When M = 1, we can derive the Rcea(z,1) and Koa(z, 1) by 
—z)? |z? — 82| 
Rce = (442)2’ Kealz, 1) = UF’ z>0. (4.24) 


There are three roots z{, = zło = 2 and zł = 8 + 6V2 of Koa(z,1) = 4, respectively. zj, = z{, = 2 is also 


the maximum value point of Kaq(z,1). Moreover, Kcaq(z,1) is an increasing function for z > 8 + 6v2. In 


conclusion, we may obtain 


1 * 
z € (0, 27] 
Zmax — 82max 


max Kcee(z,0)< 
zE[0, zmax] 


(4.25) 


ae eer 
When M > 2, we can infer from Figure 4.2 that there exists only one root depends on M z*(M) > 27 of 
Koa(z, M) = § and Keg(z, M) is an increasing function for z > z*(M), then we have 


1 , 
max Koq(z,M) = ¢ 3 e (4.26) 


2€[0,2max] Koa (zmax, M) ZE (2* (M), Zmax]- 


Therefore, we can derive that Theorem 4.2 holds. m 
We can infer from Theorem 4.2 that M has a significant impact on the convergence rates because the 


lower bound MžŽin 


increases aS Zmax = ATAmax grows. Using while loops in Matlab, we can derive M*,,, for 


a given Zmax Since M*. is a natural number. 


min 


60 T T T T 60 


1 1 
0 2000 4000 


Zmax = 


6000 8000 10000 10? 107 10° 10! 10? 10° 104 
AT max 


Figure 4.3: MŠ 
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Remark 4.1 (Comparison to other Parareal algorithms) Researchers have proved that the following 
proposition for the Parareal-Euler [380] and Parareal-TR/BDF2 [35] algorithms utilizing the backward Euler 
and TR/BDF2 (i.e. ode23 solver for ODEs in MATLAB), respectively, as the F propagators in the earlier 
papers. 


max Kuler, TR/BDF2(2; J) < V Zmax > 0, (4.27) 


1 
z€[0,zmax] 37 
which implies that for all zmax > 0, the convergence rates of Parareal algorithms p(J) < z. 

Unfortunately, not all methods yield such a uniform result. The following conclusion holds for the 
Parareal-TR and Parareal-Gauss4 algorithms [42] when utilizing the trapezoidal rule and fourth-order Gauss 
Runge-Kutta method as F-propagator. 

, if J is even and J > J; (4.28) 


min? 


max Kor, Gauasa(Z; J) < 
2€[0,2max] 


w| = 


It indicates that if the mesh ratio J is large enough, the two parareal algorithms will converge as fast as 
the Parareal-Euler method. Our Parareal-CG method behaves similarly as Parareal-TR and Parareal-Gauss4 
algorithms. 

For the explicit F-propagator forward Euler and fourth-order explicit Runge-Kutta algorithms used in 
the Parareal-fEuler and Parareal-4ERK algorithms, respectively, the convergence factor K begins to suddenly 
trend to infinity at some Zmax. Therefore, these parareal algorithms can be used to solve a very limited 
number of equations. 


KiEuler, 4ERK(Z, J) =o, SS 2 (4.29) 


The convergence factors of the three kinds of parareal algorithms are shown in Figure 4.4, 


5 Numerical Experiment 


In this section we verify the convergence and efficiency of the Parareal-CG algorithm. In all the computations, 
the initial iteration for the parareal algorithm is chosen randomly and the iteration process stops when the 
following tolerance is obtained: 


max ||uk+" — uk ||, < 107). (5.1) 
n 


Moreover, the absolute error and iteration error in all experiments are defined by 


k+1 
5 = 


k+1 


Absolute error : max ||u u(Tn)|loo. Iteration error : max ||uk*? — uë ||... (5.2) 
n n 


Example 5.1 (Near Earth satellite motion integration problem) We take a two-body motion in low 


Earth orbit (LEO) with just mutual gravitational attraction for the first example. The system of the second 


14 


202304.00993v1 


chinaXiv 


Parareal-Euler J=4 

Parareal-Euler J=20 
Parareal-TR/BDF2 J=4 | | 
Parareal-TR/BDF2 J=20 


Parareal-TR J=4 

0.9 Parareal-TR J=20 
zhsss Parareal-Gauss4 J=4 
Lit ated Parareal-Gauss4 J=20 


Parareal-forward Euler J=4 
o9 U Parareal-forward Euler J=20 
peste Parareal-4ERK J=4 

panes Parareal-4ERK J=20 


01i F 


10° 10° 10? 104 10 10? 107 10° 10' 10? 
Zmax(= ATAmax) Zmax(= AT Amax) 


Figure 4.4: The contraction factors K as a functions of Zmax, for different M. Top: Parareal- 
Euler (F=backward Euler) and Parareal-TR/BDF2 (F=TR/BDF2). Bottom left side: Parareal-TR 
(F=trapezoidal rule) and Parareal-Gauss4 (F=forth-order Gauss). Bottom right side: Parareal-fEuler 
(F=forward Euler) and Parareal-4ERK (F=forth-order explicit Runge-Kutta). 


order three-dimensional equations is 


a(t) = -gry t € [0,50], 
y"(t) = “eat t € [0,50], (5.3) 
z(t) = -agy t € [0,50], 


where x, y and z are the three coordinates in some Earth-centered inertial reference frame; r is the distance 
between the two bodies defined by 


r(t,z,y, z) = y z(t)? + y(t)? + z(t); 


u = 3.986 x 105km?/s? is the Earth gravitational constant. We formulate the equations as a first-order 
system of siz-dimensional differential equations, We formulate the equations as a first-order system of siz- 


dimensional differential equations, where the solution is uniquely determined by the initial position and 
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T T 
—©— Parareal-Euler 
—®— Parareal-TR 


—#— Parareal-Gauss4 
—A— Parareal-CG 


T T 
—©— Parareal-Euler 
—®— Parareal-TR 


—#— Parareal-Gauss4 
—A— Parareal-CG 


Absolute error 
Iteration error 


i f 1 i 1 1 1 i 1 1 f 1 i 1 
5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 
Iteration:k Iteration:k 


Figure 5.5: Comparisons of the absolute errors (left) and iteration errors (right) for the Parareal-Euler, 
Parareal-TR and Parareal-Gauss4 and Parareal-CG algorithms. 


velocity 


[x(0), y(0), z(0)] = [464.856, 6667.880, 574.231]km, 
[x’(0), y'(0), 2’(0)] = [—2.8381188, —0.7871898, 7.0830275]km/s. 


We obtain the exact solution of the equation (5.3) by solving the two-body Keplerian motion using the F and 
G approach in [33]. 


We compare the accuracy of Parareal-CG, Parareal-Euler, Parareal-TR/BDF2 and Parareal-Gauss4 al- 
gorithms revisited in Remark 4.1 for solving the example. In each coarse grid, we fix AT = 0.25 with M = 6 
for Parareal-CG algorithm and J = 6 for other parareal algorithms. Then we perform the absolute errors and 
iteration errors of the parareal algorithms in Figure 5.5. In this way, we can draw the following conclusions 


by the convergence behaviors in the figures. 


e In comparison to the Parareal-Euler and Parareal-TR, the Parareal-CG algorithm utilizes less iter- 
ations. While the Parareal-Gauss4 algorithm uses the same number of iterations as Parareal-CG 


algorithm. 
e Using the same number of points M = J = 6 in every coarse grid, the Parareal-CG algorithm has a 


higher accuracy than the three others. 


Example 5.2 (Burgers’ equation) Consider the 1D Burgers’ equation with initial and boundary condi- 


tion as following 


a aa + Mae = 0, (x,t) € [0,2) x [0,4], 

_ 2vrsin(r2) 
w(t, 0) = Foma) TELD (5.4) 
u(0,t) = u(2,t) = 0, t € [0,4], 
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in our computations, we choose a = 2 and test v = 0.05, 0.005. The exact solution of the equation is 


_ 2vmexp(—77vt) sin(rx) 


u(x,t) = (5.5) 


a + exp(—72vt) cos(mx) 


We divide the spatial domain x € [0,2) into an Ny mesh uniformly with Ax = , in this way, we have 
x; = jAx, j =0,1,---,N, — 1. Applying the fourth-order compact finite difference scheme to approximate 


gu and ot we have 
10 20 10 1 
a ae ejat) + 5 (Ent) + oq (ast) = z (ulaj42,t) — ulz, t)), 
6 Ox 3 Ox 6 Ox 2Ax (5.6) 
1 Pu 5 Ou 1 Au 1 : 
T gga C12) + 6 Daz Zitt) ED ax (zj,t) = Ar (u(xj+2,t) — 2u(x;41,t) + u(z;,t)). 
Define the vectors 
. Ou(ao,t) ulz, t) Ou(an,-1,t) i 
u= , gy ’ 
ot Ot ot (5.7) 
u = [u(xo,t), u(x1,t),--- ,u(ay,—1,t)]" 


Then combining the compact finite difference scheme with the period boundary condition, the Burgers’ 


equation (5.4) has the following spatial semidiscrete scheme 
ù + Au + u- (Agu) = 0, (5.8) 


where the coefficient matrices A; and Ag are 


-1 


a A| |-2 1 1 
b 2? = 1 22 4 
A= ->a , 
ben 1-2 1 
+ a 2 1 i <Q 
=] 
E a fea a 
23 4 =O a 
1 
ASS Sa 
3 3 3 _ ae 
3 2 3} [2 “10 


The three figures in Fig 5.6 show the dependence of the convergence rate of the Parareal-CG algorithm on 


AT, Ax, M, respectively. We can draw the following conclusions by the figure. 
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i 


Iteration Number 
sc 


a 29 Nn o © 


Iteration Number 
Iteration Number 


10" 10° 10! 102 
Ax M 


Figure 5.6: Top: The semi-log plot of the dependency of the Parareal-CG algorithm’s convergence rate on 
AT, Av = 1/4, M = 4, AT ranges from 2~3 to 278. Bottom left side: The semi-log plot of the dependency 
of the Parareal-CG algorithm’s convergence rate on Az, AT = 1/64, M = 4, Az ranges from 271 to 275. 
Bottom right side: The semi-log plot of the dependency of the Parareal-CG algorithm’s convergence rate on 
M, AT = 1/32, Ax = 1/4, M ranges from 2! to 2°. 


e Each of the v possesses robust convergence rate with respect to the change of AT. The convergence 
factor decreases as AT increases. 


e Each of the v possesses robust convergence rate with respect to the change of Ax. The convergence 
factor increases as Ax varies from small to large since the eigenvalues of A; and A» increases as Ax 


reduces. 


e The convergence rate is insensitive to the choice of M with AT = 1/32, which implies the high accuracy 


of the Chebyshev-Gauss spectral collocation method. 


e All the experiments v = 0.005 needs less iteration number than v = 0.05. 


6 Conclusion 


In this paper, we propose the Parareal-CG method for solving time-dependent differential equations, where 
the coarse propagator G is fixed by backward Euler method and the fine propagator F is chosen to be 
Chebyshev-Gauss spectral collocation method. The algorithm does have a convergence factor around 0.333, 


although the number of Chebyshev-Gauss points M needs to be somewhat large. The spectral radius of the 
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matrix A and the AT provide the lower bound of M that guarantees such a expective convergence factor. 


Some numerical experiments illustrate the accuracy and convergency of the presented algorithm. 
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